
    source("1.2.1 alpha_est.R")

    Reg = function(x.delta, x.delta.d, x.op.y, x.op.d, z, d, y, weights = NULL){
        
        if(is.null(weights)) weights = rep(1,length(y))
        
        d.rd.model = brm(d,z,x.delta.d,x.op.d,param = "RD", est.method = "MLE",
                         weights = weights)
        p.beta = ncol(x.delta.d)
        p.eta  = ncol(x.op.y)
        beta.est = d.rd.model$point.est[1:p.beta]
        eta.est  = d.rd.model$point.est[(p.beta+1):(p.beta+p.eta)]
        atanh.delta.d.est = x.delta.d %*% beta.est
        delta.d.est = tanh(atanh.delta.d.est)
        lopop.d.est = x.op.d %*% eta.est
        p0.d.est    = mapply(getProbScalarRD, atanh.delta.d.est, lopop.d.est)[1,]
        
        pa = dim(x.delta)[2];  alpha.start = rep(0,pa)
        pb = dim(x.op.y)[2];   zeta.start  = rep(0,pb)
        alpha.sol = max.likelihood(y, z, x.delta, x.op.y, alpha.start, zeta.start, 
                               weights, max.step=1000, thres=1e-6, pa, pb,
                               delta.d.est = delta.d.est)
        alpha.est = alpha.sol$par[1:pa]
        zeta.est  = alpha.sol$par[(pa+1):(pa+pb)]
        delta.est = tanh(x.delta %*% alpha.est)
        delta.y.est = delta.est * delta.d.est
        lopop.y.est = x.op.y %*% zeta.est
        p0.y.est    = mapply(getProbScalarRD, atan(delta.y.est), lopop.y.est)[1,]
        
        Delta.est = sum(weights * delta.est) / sum(weights)
        
        return(list(p0.d.est=p0.d.est, p0.y.est=p0.y.est, delta.d.est=delta.d.est,
                    delta.est=delta.est, Delta.est=Delta.est, alpha.est=alpha.est,
                    beta.est=beta.est, zeta.est=zeta.est, eta.est=eta.est))
        
    }